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THE RADIAL DISTRIBUTION OF GALAXIES IN ACDM CLUSTERS 
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ABSTRACT 

We study the radial distribution of subhalos and galaxies using high-resolution cosmological simulations 
of galaxy clusters formed in the concordance ACDM cosmology. In agreement with previous studies, 
we find that the radial distribution of subhalos is significantly less concentrated than that of the dark 
matter, when subhalos are selected using their present-day gravitationally bound mass. We show that the 
difference in the radial distribution is not a numerical artifact and is due to tidal stripping. The subhalos 
in the cluster core lose more than 70% of their initial mass since accretion, while the average tidal mass 
loss for halos near the virial radius is ~ 30%. This introduces a radial bias in the spatial distribution 
of subhalos when the subhalos are selected using their tidally truncated mass. We demonstrate that the 
radial bias disappears almost entirely if subhalos are selected using their mass or circular velocity at the 
accretion epoch. The comparisons of the results of dissipationless simulations to the observed distribution 
of galaxies in clusters arc therefore sensitive to the selection criteria used to select subhalo samples. Using 
the simulations that include cooling and starformation, we show that the radial distribution of subhalos is 
in reasonable agreement with the observed radial distribution of galaxies in clusters over the entire radial 
range probed by the simulations, 0.1 < i?/i?200 < 2.0, if subhalos are selected using the stellar mass of 
galaxies they harbor. The radial bias is minimized in this case because the stars are located in the centers 
of dark matter subhalos and are tightly bound. The stellar mass of an object is therefore approximately 
conserved as the dark matter is stripped from the outer regions. Nevertheless, the concentration of the 
radial distribution of galaxies is systematically lower than that of the dark matter. Finally, we find that 
the baryon dissipation does enhance the survival of subhalos, especially in the cluster cores. However, 
the effect is relatively weak because the baryon dissipation affects the density distribution only at small 
radii. 

Subject headings: cosmology: theory-galaxies: formation- methods: numerical 



1. INTRODUCTION 

Understanding the spatial distribution of galaxies and 
matter is one of the main goals of observational and theo- 
retical cosmology. At large scales, the galaxy and matter 
distributions on average are thought to be simply related 
by a constant scale-independent bias, with the constant 
depending on the galaxy mass (e.g., Mo & White 1996; 
Scherrer & Weinberg 1998). Clustering of galaxies at small 
scales (r < 3/jT 1 Mpc), on the other hand, is thought to 
be determined by the abundance and radial distribution 
of galaxies within dark matter halos (e.g., see Cooray & 
Sheth 2002, for a recent review). 

While the abundance of galaxies in halos of different 
mass have recently been extensively quantified in cosmo- 
logical simulations (e.g., Seljak 2000; Guzik & Seljak 2001; 
White et al. 2001; Berlind et al. 2003; Kravtsov et al. 
2004a), their radial distribution is less well understood. 
Theoretical analyses often simply assume that galaxies 
trace the radial profile of dark matter. However, it is easy 
to imagine that dynamical friction, tidal disruption, and 
morphological transformation could modify the distribu- 
tion of galaxies with respect to matter. In addition, the 
radial distribution of galaxies should clearly be sensitive to 
how galaxies are selected. For example, the distribution of 
early type red galaxies in clusters is more centrally concen- 
trated than that of late type blue galaxies (e.g., Goto et al. 
2004) . New cluster samples based on large galaxy surveys 
(Carlberg et al. 1997; van der Marel et al. 2000; Lin et al. 
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2004; Miller et al. 2004) should shed light on these issues, 
but the interpretation of observational results will require 
comparisons with theoretical expectations. 

Properties of galaxy-size dark matter halos in groups 
and clusters (or the subhalos) have been the subject of 
many recent studies, which used a new generation of high- 
resolution dissipationless simulations not affected by the 
"overmerging" problem (Ghigna et al. 1998; Tormen et al. 
1998; Klypin et al. 1999; Coh'n et al. 1999; Okamoto & 
Habe 1999; Colin et al. 2000; Ghigna et al. 2000; Springel 
et al. 2001; Taffoni et al. 2003; De Lucia et al. 2004; Tor- 
men et al. 2004; Kravtsov et al. 2004b; Diemand et al. 
2004; Gao et al. 2004b; Reed et al. 2004). One of the main 
results is that the radial distribution of subhalos is less 
concentrated than that of dark matter (Ghigna et al. 1998; 
Colin et al. 1999; Ghigna et al. 2000; Springel et al. 2001; 
De Lucia et al. 2004; Gao et al. 2004a). The subhalos also 
appear to have a significantly more extended and shallower 
radial distribution compared to the observed distribution 
of galaxies in clusters, a point recently emphasized by Die- 
mand et al. (2004) and Gao et al. (2004a). 

Although questions were raised about remaining over- 
merging in the inner regions of halos (Taylor et al. 2003), 
Diemand et al. (2004) used convergence tests to show that 
the subhalo distribution is not affected by resolution but is 
due to real physical merging of subhalos in dissipationless 
simulations (see also § 4 below). They also argued that the 
physical overmerging in dissipationless simulations results 
in a factor of two underestimate of the number of galactic 
subhalos in cluster-size systems. 
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One of the obvious omissions in dark matter only simu- 
lations is the lack of dissipation and starformation. After 
all, the observed samples of galaxies are most often se- 
lected using their luminosity or an observational proxy for 
stellar mass, not a total gravitationally bound mass. In 
addition, when baryons cool and condense in the cores of 
dark matter halos, they increase the inner density of their 
host, which can make it more resistant to tidal disrup- 
tion (White & Rees 1978). It is therefore important to 
test whether inclusion of cooling and starformation, the 
processes critical for realistic modeling of galaxy forma- 
tion, results in significantly different galaxy populations 
in clusters. 

Although a number of studies during the last decade 
used gasdynamics cosmological simulations with cooling 
to study galaxy clustering on small scales (e.g., Katz et al. 
1992; Pearce et al. 1999; Blanton et al. 2000; Yoshikawa 
et al. 2001; Weinberg et al. 2004), only a few analyses di- 
rectly addressed the radial distribution of galaxies in halos. 
Frenk et al. (1996) used simulations with cooling and star- 
formation to study distribution and dynamics of galax- 
ies in clusters. They found that stellar systems formed 
in the simulations were more centrally concentrated than 
dark matter. This, however, was the case only in the runs 
with starformation. When gas was allowed to cool with- 
out forming stars, the simulations appeared to suffer from 
the "overmerging" problem. More recently, Berlind et al. 
(2003) analyzed the radial distribution of galaxies, identi- 
fied as dense baryonic clumps, in halos of different mass 
in SPH simulations. These authors found that galaxies se- 
lected to have baryonic masses above a certain threshold 
have a somewhat more extended distribution compared to 
the dark matter. 

A number of studies used a hybrid approach of combin- 
ing high-resolution A-body simulations with semi-analytic 
models of galaxy formation to study galaxy populations in 
clusters and groups (Springel et al. 2001; Okamoto & Na- 
gashima 2001; Kravtsov et al. 2004b; Gao et al. 2004a). 
All of these models imply that radial distribution of halos 
is more extended than that of the luminous galaxies. By 
necessity, semi-analytic studies make certain assumptions 
about dynamical evolution of the stellar component, not 
included directly in simulations. For example, Springel 
et al. (2001) and Gao et al. (2004a) find that for a large 
fraction of objects in the cluster core the stellar systems 
may survive as distinct galaxies, even after their DM halos 
are disrupted. Kravtsov et al. (2004b) used the tidal force 
measured in simulations to estimate the amount of tidal 
heating and mass loss of galaxies semi-analytically. 

It is clearly important to address the radial distribution 
of subhalos and galaxies in self-consistent high-resolution 
cosmological simulations of cluster formation. In this pa- 
per, we present such a study. We compare results of cluster 
simulations with gasdynamics, cooling, and starformation 
to the dissipationlcss simulations started from the same 
initial conditions. We test the convergence of our results 
by comparing results of two simulations of the same cluster 
with an order of magnitude different resolution. 

The paper is organized as follows. In the following two 
sections, we describe the numerical simulations and halo 
finding algorithm used in our analysis. In § 4 we present 
the convergence tests. We present the results of dissipa- 



tionlcss simulations in § 5 and gasdynamics simulations in 
§ 6. We summarize our results and conclusions in § 7. 

2. SIMULATIONS 

In this study, we analyze high-resolution cosmological 
simulations of eight group and cluster-size systems in the 
"concordance" flat ACDM model: Q m = 1 - A = 0.3, 
n h = 0.043, ft = 0.7 and cr 8 = 0.9, where the Hubble 
constant is defined as 100ft- km s _1 Mpc -1 , and as is the 
power spectrum normalization on 8ft" 1 Mpc scale. The 
virial masses of cluster systems we consider range from 
« 7 x lO^ft" 1 M Q to 3 x lO^ft" 1 Mq. The simulations 
were done with the Adaptive Refinement Tree (ART) A- 
body+gasdynamics code (Kravtsov 1999; Kravtsov, Klypin, 
& Hoffman 2002), an Eulerian code, which uses the adap- 
tive refinement in space and time, and (non-adaptive) re- 
finement in mass (Klypin et al. 2001) to reach the high 
dynamic range required to resolve cores of halos formed in 
self-consistent cosmological simulations. 

For each cluster, we analyze two sets of simulations 
started from the same initial conditions but with different 
physical processes included. The first set of simulations 
followed dynamics of dark matter only. To study effects of 
resolution, one of the clusters was re-simulated with eight 
times more particles and with higher spatial resolution 
(Tasitsiomi et al. 2004). We will be denoting lower and 
higher resolution simulations of this cluster as LR and HR 
runs. Both LR and HR runs used a 256 3 uniform root grid 
covering the computational box of 80ft. -1 Mpc, but differ- 
ent numbers of particles and refinements. The LR run has 
an effective mass resolution of Mp R = 3.16 x 10 8 /i _1 M©, 
corresponding to 512 3 particles in the box, and reached 8 
levels of refinement, corresponding to the size of the high- 
est refinement level cell of 1.2ft, -1 kpc. The HR run has 
eight times more particles (m™ = 3.95 x 10 7 ft _1 Mq). In 
this simulation, the smallest cell size reached was 0.6ft" 1 
comoving kpc. 

The gasdynamics simulations started from the same ini- 
tial conditions as the A-body runs used a 128 3 uniform 
grid and 8 levels of mesh refinement, which corresponds 
to the dynamic range of 128 x 2 8 = 32768 and peak res- 
olution of 80/32,768 « 2.44ft" 1 kpc. Only the region of 
~ 10ft" 1 Mpc around the cluster was adaptively refined, 
the rest of the volume was followed on the uniform 128 3 
grid. The dark matter particle mass in the region around 
the cluster was 2.7 x 10 8 ft, _1 M©, while other regions were 
simulated with lower mass resolution. Note that although 
the mass resolution of the gasdynamics and the LR N- 
body run is similar, the spatial resolution in the latter is 
higher. In addition, more aggressive refinement criteria 
were used in the A-body runs compared to the gasdy- 
namics simulations. Some differences in the completeness 
of subhalo samples are therefore expected. 

The gasdynamics simulations included dynamics of gas 
and collisionless DM and several physical processes criti- 
cal to various aspects of galaxy formation: star formation, 
metal enrichment and thermal feedback due to the super- 
novae type II and type la (SNII/Ia), self-consistent advec- 
tion of metals, metallicity- and density-dependent cooling 
and UV heating due to cosmological ionizing background 
using cooling and heating rates tabulated for the temper- 
ature range 10 2 < T < 10 9 K and a grid of densities, 
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metallicities, and UV intensities using the Cloudy code 
(ver. 96b4, Ferland et al. 1998). The cooling and heat- 
ing rates take into account Compton heating/cooling of 
plasma, UV heating, atomic and molecular cooling. Star 
formation in the cluster simulations is implemented us- 
ing the observationally- motivated recipe (e.g., Kennicutt 
1998): p. = Pgas/i*, with i, = 4x 10 9 yrs. Stars are al- 
lowed to form in regions with temperature T < 2 x 10 4 K 
and gas density n > 0.1 cm -3 . A more detailed description 
of these simulations will be presented elsewhere. 

3. HALO IDENTIFICATION 

A variant of the Bound Density Maxima halo finding 
algorithm (Klypin et al. 1999) is used to identify halos 
and the subhalos within them. The details of the algo- 
rithm and parameters used in the halo finder can be found 
in Kravtsov et al. (2004a). The main steps of the algo- 
rithm are identification of local density peaks (potential 
halo centers) and analysis of the density distribution and 
velocities of the surrounding particles to test whether a 
given peak corresponds to a gravitationally bound clump. 
More specifically, we construct density, circular velocity, 
and velocity dispersion profiles around each center and it- 
eratively remove unbound particles using the procedure 
outlined in Klypin et al. (1999). We then construct final 
profiles using only bound particles and use these profiles 
to calculate properties of halos, such as the circular ve- 
locity profile Vd r c(r) = y/GM(< r)/r and compute the 
maximum circular velocity Vm ax . 

If the center of a halo does not lie within a larger system, 
we consider the halo to be an isolated or host halo. In this 
case we use the virial radius, r v j r , defined as the radius en- 
closing an overdensity of 180 with respect to the mean den- 
sity of the Universe 2 . For halos located within the virial 
radius of a larger host halo (the subhalos), we define the 
outer boundary at the truncation radius, r t , at which the 
logarithmic slope of the density profile constructed from 
the bound particles becomes larger than —0.5 as we do not 
expect the density profile of the CDM halos to be flatter 
than this slope. Throughout this paper, we will denote the 
minimum of the virial mass and mass within r t , simply as 
M. In the gasdynamics simulations, the total mass of an 
object includes the mass of dark matter, gas, and stars. 
For each system we also estimate the stellar mass, M*, 
within the truncation radius. 

4. CONVERGENCE TESTS 

Before we proceed with our analysis, it is important to 
test the effects of resolution and determine the smallest 
mass of objects that are resolved and identified reliably in 
our simulations. Figure 1 shows comparison of the cumu- 
lative mass functions of subhalos in the two A-body and 
gasdynamics simulations of the same cluster. We present 
separate comparisons for the objects within the virial ra- 
dius and within the inner half of the virial radius to show 
any dependence on radius. 

Note that small differences both between gasdynam- 
ics and A-body simulations and between low- and high- 

2 We will use this definition of the virial radius throughout the paper. 
In § 6.3, to compare the radial distributions of galaxies in simulations 
and observations we will use -R200 > defined as the radius enclosing the 
overdensity of 200 with respect to the critical density of the Universe. 
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Fig. 1. — Cumulative mass function of subhalos enclosed within 
r v ; r (top) and 0.5r v ; r (bottom) of the galaxy cluster with M v ; r = 
2.4 X 10 14 h — 1 Mq. Three lines indicate gasdynamics run with cool- 
ing and starformation (solid), LR (dotted), and HR (dashed) A-body 
runs. The mass functions of all three runs converge for A p > 80, 
which corresponds to M/M V \ T > 10~ 4 for hydro and LR runs and 
M/M viT > 10~ 5 for the HR A-body run. 



resolution dissipationless simulations are expected even at 
masses where convergence is reached. The former can arise 
due to different time integration in the A-body and hydro 
codes (see, Frenk et al. 1999), while the latter can be due 
to the fact that the HR run includes small-scale modes 
absent in the LR run. With this in mind, comparisons of 
the mass functions in the LR A-body run and gasdynam- 
ics run agree well with the results of the HR A-body run 
for A p > 80 or M/M v ; r > 1 x 10~ 4 . Similar comparison 
of the velocity functions of subhalos in the LR A-body 
run agree well with the results of the HR A-body run for 
Vm/Kir > 0.1 or V m > 100 km/s for the cluster halo with 
Vyit — 1065 km/s. We will therefore only consider objects 
with masses and maximum circular velocities above these 
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Fig. 2. — Convergence test for the radial distribution of subha- 
los in the same cluster as in Figure 1. Number density profiles of 
subhalos normalized to the number density within the virial radius, 
n(r)/(n, v ; r }, in the HR (dotted) and LR (dashed) TV-body runs for 
subhalos with masses M/M v i T > 10~ 4 . The solid line shows the 
dark matter profile normalized to the virial overdensity. The figure 
shows that the convergence for the radial distribution of subhalos is 
reached, as the profiles are similar in LR and HR runs. 



limits in our analyses. We find no radial dependence of 
the slope of the cumulative mass and velocity functions. 
These results are in agreement with the convergence stud- 
ies of Diemand et al. (2004) and Gao et al. (2004b). 

Figure 2 shows the radial number density profiles of sub- 
halos with total masses > 10 -4 M V - 1T in the LR and HR 
iV-body runs. The profiles are obtained by averaging over 
five (forty) outputs between z = 0.25 and z = for the LR 
(HR) iV-body runs, respectively 3 . For comparison of dark 
matter and subhalos, we normalize the profiles to the mean 
density within the virial radius, (n„j r ). The figure shows 
that the profiles agree at all probed radii. The convergence 
is real and the subhalo profiles match without normaliza- 
tion (in other words, for the same mass threshold (n V i r ) 
is the same) . This shows that both the normalization and 
shape of the radial profile of subhalos have converged for 
the current resolution of the LR N-body run. We therefore 
conclude that the difference between matter and subhalo 
distributions is real and is not due to numerical overmerg- 
ing, in agreement with Diemand et al. (2004) and Gao 
et al. (2004b). This is inconsistent with the semi-analytic 
models of Taylor et al. (2003), which predict a consider- 
ably more radially concentrated distribution of subhalos 
than is observed in simulations. 

5. DISSIPATIONLESS SIMULATIONS 

In this section, we present the analyses of the dissipa- 
tionless simulations to study the radial distribution of sub- 
halos and how it depends on the selection criteria used to 
define subhalo sample. Figure 3 shows the radial number 

3 The difference in the number of epochs used in averaging is due to 
fewer outputs saved in the LR run. 



density profiles normalized to the mean number density of 
subhalos within the virial radius in the HR iV-body run. 
The four panels show the radial distribution of subhalos 
with different selection criteria. The top panels show the 
profile where the subhalos are selected using the present- 
day values of the subhalo mass, M , and the maximum 
circular velocity, V^. The bottom panels show the profiles 
where the values at the time of accretion onto the cluster 
are used. These values are obtained from the evolution- 
ary tracks of each subhalo constructed as described by 
Kravtsov et al. (2004b) . We define the accretion epoch as 
the time when a halo's distance to the most massive clus- 
ter progenitor first becomes smaller than 1.5 of its virial 
radius (at this epoch). 

The profiles are constructed by stacking forty outputs 
between z = 0.25 and z — to increase the statistics 
and reduce fluctuations. Note that this cluster does not 
evolve significantly between z = 0.25 and z — 0. We also 
checked that the results are not biased by the stacking 
procedure. Throughout the paper, the error bars indicate 
la Poisson error defined by the number of objects in each 
bin. Here and throughout, we normalize all profiles to the 
mean density within the virial radius because this allows 
unambiguous comparison of radial distributions of dark 
matter and of halos of different masses. We do not include 
the central galaxy in the first bin of the profiles. 

We show each selection using three values of threshold 
mass or circular velocity. The values of the circular ve- 
locity thresholds are chosen to match the corresponding 
mass thresholds approximately. The shape of the profiles 
does not depend on the subhalo mass and is significantly 
shallower than the dark matter distribution for all mass 
thresholds. This is consistent with the results of Diemand 
et al. (2004). 

The more extended distribution of subhalos compared to 
DM in this case is simply due to the fact that the subhalos 
in the inner regions have on average suffered larger tidal 
mass loss than the halos near the virial radius. Figure 4 
shows the fractional mass loss and the change in V^ ax ex- 
perienced by each subhalo since the epoch when it reached 
the maximum mass, M max , in its evolution. The scatter 
for individual objects is large and non-Gaussian and is 
due to the wide distribution of accretion times and orbital 
parameters of subhalos. Nevertheless, there is a clear av- 
erage trend of increasing mass loss at smaller radii. For 
instance, the halos within 0.3r v i r on average have lost more 
than 70% of their mass since accretion. At r > 0.5r v ; r the 
halos on average lose only < 40% of their mass. Note 
also that the tidal mass loss is also accompanied by a slow 
decrease in Vmax (Kravtsov et al. 2004b, see also Tormen 
et al. 2004). Mass-based selection thus biases the sub- 
halo sample to large radii where the tidal stripping have 
depleted the number of objects in a given mass range to a 
lesser degree. 

Figure 3 shows that the bias introduced by selecting 
subhalo samples using circular velocity thresholds, V^, is 
smaller than in the selection based on mass. The bias is ex- 
pected in this case because circular velocity evolves when 
a halo experiences tidal loss as V m oc M°- 2 ~ - 25 (Kravtsov 
et al. 2004b). The average relation between circular ve- 
locity and mass for isolated halos is V m oc M - 3 (e.g., 
Avila-Reese et al. 1999; Bullock et al. 2001). The mass 
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Fig. 3. — Radial distribution of subhalos for samples with different selection criteria. Top panels: the subhalos are selected using the 
present-day values of the subhalo total mass, M°, (left) and the maximum circular velocity, (right). Bottom panels: the subhalos are 
selected using the subhalo mass, M acc , (left) and the maximum circular velocity, V^ cc (right) at the time of accretion. In each panel the 
subhalos in the three samples were selected using either the minimum mass thresholds of M/M v i T > 10~ 4 (dotted), > 10~ 3 (dashed) and 
> 2 X 10 — 3 (long-dashed) or corresponding values of the circular velocity: Vm > 60kms — 1 (dotted), > 128 km s -1 (dashed) and > 162 km s -1 
(long- dashed). The profiles are constructed by stacking 40 outputs of the HR TV-body run between z = 0.25 and z = 0. The error bars are 
the lcr Poisson errors, and the central galaxy is not included in the first bin of the profiles. 



evolution would thus tend to shift the halos off the rela- 
tion. Indeed we find that the normalization of the M — V m 
relation changes for the smaller cluster-centric distances, 
although the slope of the relation is approximately con- 
stant. The slower evolution of the maximum circular ve- 
locity compared to mass means that when it is used to 
select subhalos the radial bias is smaller. 

If the radial bias is due solely to the varying amount of 
mass loss at different radii, we can expect that it should 
be considerably smaller if we select subhalos using their 
mass or circular velocity at the time of accretion (i.e., not 
modified by the tidal stripping yet). The bottom panels of 
Figure 3 show the radial distribution of subhalos for such 
selections. It is clear that the subhalo radial profiles in 
this case are rather similar to the dark matter profile at 
r > 0.2r v j r . The mass measured at the accretion, thus, 
provides a "label" for each object, which is not affected by 



the subsequent tidal mass loss and evolution. Note that we 
can still expect a certain bias because a larger fraction of 
subhalos is fully disrupted by tides at small radii. We can- 
not compensate for the absence of the disrupted subhalos 
with the change of selection criteria for halos which sur- 
vived to z = 0. We checked that the results are insensitive 
to the definition of time of accretion, as long as subha- 
los are selected based on their mass and V m well before 
entering into the cluster virial region. 

Our results therefore indicate that the bias in the radial 
distribution of subhalos with respect to dark matter arises 
simply due to the radial dependence of the effects of tidal 
stripping. It is minimized if we use a property of halos 
unaffected by the evolution. In dissipationless simulations 
such properties can only be computed at the actual time 
of accretion if the evolution of a halo is traced in time. 
Unfortunately, the halo properties at the time of accretion 




Fig. 4. — The fractional mass loss and the change in V,^ ax be- 
tween the epochs when each subhalo has reached the maximum 
mass, M max , and z = as a function of the distance to the cluster 
center. The dots show the mass loss experienced by individual sub- 
halos, while the solid circles show the median values in logarithmic 
radial bins. The triangles shows the same for the change in the fifth 
power of Vmax- Halos at smaller radii have on average experienced 
considerably more stripping than subhalos at larger radii. The tidal 
mass loss is also accompanied by the decrease in V ma x- 

are not observable in reality. In the next section, we show 
that the stellar mass is an observable with similar prop- 
erties. The stellar mass of an object should change little 
until it is almost fully destroyed by tides because stars are 
located in the centers of halos and are tightly bound. We 
can expect, therefore, that stellar mass should correlate 
well with the mass and circular velocity at the time of ac- 
cretion and that the distribution of galaxies selected using 
stellar mass (or luminosity) should be similar to that of 
dark matter. 

6. GASDYNAMICS SIMULATIONS 

6.1. Radial distribution of subhalos 

In this section we consider the effects of gas cooling and 
starformation on the abundance and radial distribution 
of subhalos. Figure 5 compares the radial distribution of 
subhalos with M/M vit > 1CT 4 and M/M vir > 2 x 10 -3 in 
the A^-body and gasdynamics simulations of the same clus- 
ter. These profiles are constructed by stacking five (nine) 
outputs between z = 0.25 and z = for TV-body (gas- 
dynamics) simulation. The radial distribution of subhalos 
in the iV-body and gasdynamics simulations are remark- 
ably similar. 4 This indicates that the gas cooling has rel- 

4 Although we plot profiles normalized to (n v ir) to compare to the 
DM profile, the values of (n v ; r ) are the same in the gasdynamics and 
dissipationless run, as can be seen in Figure 1. The unnormalized 
profiles do match in the same way as in Figure 5. 
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Fig. 5. — Radial profiles n(r) / (n v { r ) in gasdynamics (dotted) and 
A^-body (dashed) runs for halos with M/M v j r > 10 -4 (top) and 
M/M v ; r > 2 X 10~ 3 (bottom). The solid line shows the dark matter 
profile in the gasdynamics run. The radial distribution of subhalos 
is remarkably similar for these two runs, regardless of the mass cut. 
This indicates that gas cooling in our simulations does not signifi- 
cantly affect the survival of subhalos and their radial distribution. 

atively little effect on the survival and radial distribution 
of subhalos. Although cooling can clearly help survival 
of subhalos in lower resolution simulations, the halos in 
our A'-body simulations are apparently sufficiently dense 
to avoid premature disruption. Given the results of the 
recent studies of the causes of the overmerging problem 
(e.g., Moore et al. 1996; Klypin et al. 1999), this is not 
surprising. Gas cooling significantly affects mass distri- 
bution only in the inner few percent of the virial radius, 
while the survival is largely determined by the halo density 
within the radius of V max (« 2r s ), just outside the affected 
regions. We should note that survival in the cluster core 
is probably enhanced by cooling. For the subhalos above 
the mass threshold of MjM w \ T > 10~ 4 , the radial profile in 
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Fig. 6. — The mass of stars vs. the total mass of subhalos lo- 
cated within r v i r (all circles) and 0.3r v j r (filled circles) of the host 
cluster. The plot shows that M* scales with Mtot with large scat- 
ter, particularly for subhalos with Mtot < 10 1 h~ Mq. Note that 
the normalization of the relation changes in the inner regions of the 
cluster, probably due to the higher average tidal mass loss of the 
objects located there. 
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Fig. 7. — Radial profiles n(r)/(n v i v ) of subhalos selected using 
the stellar mass thresholds of M„ > 10° (dotted), 10 9 (dashed) and 
10 10 h~ 1 Mq (long-dashed) in the gasdynamics run. Note that the 
radial distribution of subhalos with the M„ -selection does not de- 
pend strongly on stellar mass threshold and is close to the radial 
density profile of dark matter. 



the gasdynamics simulation is steeper than in the dissipa- 
tionless run at r < 0.2r v j r . This, however, does not change 
the overall difference between subhalo and DM profiles at 
larger radii. 

6.2. Radial distribution of galaxies 

In reality, we observe galaxies, not subhalos, and it 
would be interesting to study the distribution of galaxies in 
simulations directly. To interpret the results, we first need 
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Fig. 8. — Radial profiles n(r)/{n v i T ) of subhalos in the gasdy- 
namics run selected using the stellar mass (dotted) and correspond- 
ing total mass (dashed). The solid line shows the radial density 
profile of dark matter in the gasdynamics simulation. We show 
the profiles for two different mass thresholds: M* > 10 9 h~ 1 MQ 
and Mtot > 3 X 10 10 h' 1 M Q (top panel), M* > 10 10 h- 1 M e and 
Mtot > 2 X 10 11 /i~ 1 Mq (bottom panel). Selection using stellar mass 
results in a significantly steeper radial profile. 



to understand how stellar mass relates to the total mass 
of the host subhalo. We compute the stellar mass, M*, 
as the mass enclosed within the tidal truncation radius of 
each subhalo. Figure 6 shows that the stellar mass of a 
galaxy, M*, correlates with the subhalo mass, M* cx M" ot . 
The slope of the correlation ranges in a ~ 1 — 1.5. How- 
ever, there is significant scatter, which becomes increas- 
ingly larger for the smaller subhalo masses. Note that the 
galaxies in the core of the cluster (r < 0.3r V i r ) populate 
the region of the plot near the upper envelope of points. 
In other words, for a given stellar mass, the total subhalo 
mass on average becomes smaller in the inner regions of 
cluster. 
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Fig. 9. — Radial profiles of subhalos n(r) / (n v lr) f° r the eight simu- 
lated clusters at z = 0. Here, subhalos are selected using their stellar 
(top) and corresponding total (bottom) masses of M, > 10 10 h~ 1 Mq 
and Mtot > 2 X 10 11 h. -1 Mq respectively. The solid line is the den- 
sity profile of dark matter averaged over all clusters. The figure 
shows that the distribution of subhalos selected using M* follows 
the dark matter profile at 0.1 < r/r v ; r < 1 reasonably well. How- 
ever, there is significant variation in the inner slope of the galaxy 
number density profiles among different clusters; for two out of the 
eight clusters the profile in the center is considerably shallower than 
that of the DM. 

This systematic trend likely arises because the subhalos 
preferentially lose DM mass from their outer regions, while 
the tidal mass loss of more tightly bound stars is consider- 
ably smaller. In this CclSC, clS the subhalos are stripped and 
total mass decreases, the stellar mass stays approximately 
constant resulting in the evolution of points horizontally 
to the left in the M to t — M* plane. As we show in Fig- 
ure 4, the subhalos in the cluster core on average have 
experienced larger tidal mass loss than the halos near the 
virial radius, which explains the systematic shift between 
solid and open points in Figure 6. The amount of tidal 



mass loss experienced by each object depends on its epoch 
of accretion, orbital parameters, and to some extent on 
its internal structure. The evolutionary differences among 
subhalos would explain the large scatter of the M tot — M* 
correlation. The maximum circular velocity of subhalos in 
simulations with cooling changes much more slowly than 
the tidally bound total mass. Consequently, we find that 
the correlation of Vm ax and M* is considerably tighter than 
the Mtot — M* correlation (Nagai & Kravtsov 2004, in 
preparation) . 

Figure 7 shows the radial distribution of galaxies se- 
lected using three different stellar mass cuts: M* > 10 8 , 
10 9 , and lO lo /i _1 M . The figure shows that the radial dis- 
tribution of galaxies with the stellar mass selection does 
not depend strongly on stellar mass and is close in shape 
to the dark matter profile at r > 0.1 — 0.2r V j r . Note that in 
the simulations with cooling, the DM profile at r < 0.1r V j r 
is somewhat steeper than the corresponding profile in the 
dissipationless simulation (Gnedin et al. 2004). It is clear 
that M*-based selection results in a radial distribution of 
galaxies which is considerably steeper than the radial pro- 
file of subhalo samples selected using total mass. 

Figure 8 shows this more clearly using direct comparison 
of the radial distribution of objects selected using stellar 
and corresponding total mass. The thresholds in M to t are 
chosen to correspond approximately to the M* thresholds 
using the M tot — A/* correlation (Figure 6). The figure 
shows that the radial distribution of objects in M*-selected 
samples is significantly steeper than the profiles of objects 
selected by the total mass. As expected, the difference 
is similar to the difference of profiles of subhalos selected 
using total mass at z = and at the time of accretion, 
seen in § 5. The stars are located near the centers of 
DM subhalos and are tightly bound. The stellar mass 
thus should not evolve significantly even if an object sheds 
a large fraction of its DM halo (e.g., Gnedin 2003, and 
references therein). Our results show that this is indeed 
the case. 

Since the results so far are based on the analyses of one 
cluster, we investigate if there is a variation in the dis- 
tribution of galaxies among different clusters. Figure 9 
shows the 3D radial distribution of galaxies for the sample 
of eight clusters at z = 0. We find that for all clusters the 
radial distribution of galaxies selected based on the stel- 
lar mass, M„, is steeper than that of the samples selected 
using the total mass, M tot . The figure also illustrates sig- 
nificant variation in the inner slope of the profile among 
different clusters. The radial profiles of galaxies in most 
clusters are close to the dark matter profiles at all radii, 
but there are some that show significant flattening of the 
galaxy number density profile at r < 0.3r vir . We find that 
the variation among profiles does not decrease if we use 
clusters at the epochs when they are most relaxed. This 
indicates that the flattening of the galaxy number den- 
sity profile is not related to the host halo mass or to the 
dynamical state of clusters but may be related to their 
evolutionary histories. 

It would also be interesting to ask whether the concen- 
tration of the radial profile of galaxies correlates with the 
concentration of the DM distribution. In Figure 10 we 
compare the concentration parameters of galaxies, CQai, 
and dark matter, cdm ; for the eight simulated clusters at 
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Fig. 10. — Comparison of the concentration parameters of galax- 
ies, CQ a i, and dark matter, com, for the eight simulated clusters 
at z = 0. CQ a i[2D], is obtained by fitting the NFW profile to the 
projected (2D) galaxy distribution, whereas cdm[3D], is the concen- 
tration parameter of dark matter obtained by fitting the 3D dark 
matter profile. 



z = 0. For galaxies we use the concentrations obtained by 
fitting the projected Navarro-Frenk- White profile (here- 
after NFW, Navarro et al. 1997) to the projected (2D) 
galaxy distribution, because the projected profiles are di- 
rectly observable. For the dark matter, on the other hand, 
we use the true concentration obtained by fitting the NFW 
to the 3D dark matter distribution. The details of the fit- 
ting procedure are described in the Appendix. 

The concentration parameters exhibit significant varia- 
tion among different clusters. The concentration parame- 
ters of dark matter span a wide range from 6 to 16. Note 
that the concentration of dark matter in the current simu- 
lations may be systematically larger than the prediction of 
the dissipationless simulations (Navarro et al. 1997; Bul- 
lock et al. 2001; Eke et al. 2001; Tasitsiomi et al. 2004) 
because the baryon dissipation steepens the inner profiles 
of dark matter (Gnedin et al. 2004). Two of the eight clus- 
ters have very similar best fit concentrations of galaxy and 
DM profiles. There is also one cluster for which the con- 
centration of galaxy profile is larger than that of dark mat- 
ter. Nevertheless, the concentration parameter of galaxies 
in the majority of clusters is smaller than that of dark 
matter. These results indicate that fits to the projected 
radial profiles of galaxies may be a poor probe of the DM 
concentrations. 

6.3. Comparison to observations 

Here we compare the distribution of galaxies and sub- 
halos in the ACDM clusters in our simulated sample to 
the observed galaxy distribution. Figure 11 shows the pro- 
jected radial distribution of dark matter, galaxies and sub- 
halos averaged over the eight simulated clusters at z = 0. 
The profiles are constructed by projecting spherical vol- 
ume with radius of 3r V i r centered on the cluster center 
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Fig. 11. — Projected radial distribution of galactic subhalos 
S(iJ)/S„i r averaged over the eight simulated clusters at z = and 
three orthogonal projections. We show the radial distribution of 
subhalos selected using the mass thresholds of M, > 10 w h,- 1 M Q 
{dotted) and M tot > 2 X 10 11 /i _1 A/ Q {dashed). The solid line is the 
average projected profile of dark matter in the gasdynamics simu- 
lation. The solid circles are the average radial profile of galaxies 
in clusters in the CNOC survey (Carlberg et al. 1997), and the 
open triangles are the average profile in the 2MASS cluster survey 
(Lin et al. 2004). The data points are scaled arbitrarily. Note that 
the distribution of M* -selected subhalos is consistent with the ob- 
served distribution of galaxies over the entire range of probed radii: 
0.1 < -R/-R200 < 2.0. 



and averaging over all simulated clusters and three orthog- 
onal projections of each cluster. The profiles are normal- 
ized in units of the mean surface density of dark mat- 
ter, E vir = M vir /irrl ir , and galaxies/subhalos, A vir /7rr^ ir , 
where M V i T and A v ; r are the mass and the number of ob- 
jects enclosed within the sphere with the virial radius r v i r 
respectively. The data points are from the CNOC (Carl- 
berg et al. 1997) and 2MASS (Lin et al. 2004) cluster sur- 
veys. The observed galaxy profiles are scaled arbitrar- 
ily and the 2MASS survey profile has been re-binned to 
smaller number of bins for clarity. 

The figure shows that the projected distribution of galax- 
ies in both simulations and observations is more concen- 
trated than that of subhalos selected using the total mass. 
The radial profile of galaxies is somewhat shallower than 
the average dark matter distribution. Nevertheless, the 
figure shows that the distribution of galaxies in the simula- 
tion is in reasonably good agreement with the observed dis- 
tribution of galaxies for the entire probed range of scales: 
0.1 < i?/i?200 < 2.0. The agreement with the data is con- 
siderably better if two clusters with very flat inner profiles 
(see Figure 9) are excluded. 

7. DISCUSSION AND CONCLUSIONS 

We presented a study of the radial distribution of subha- 
los and galaxies using high-resolution cosmological simula- 
tions of galaxy clusters formed in the concordance ACDM 
cosmology. In this paper, we first analyze the dissipa- 
tionless simulations to study properties and processes that 
govern the evolution of subhalos. We then analyze cluster 
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simulations with gasdynamics, cooling, and starformation 
started from the same initial conditions to study the effect 
of these processes on the spatial distribution of subhalos 
and galaxies. 

In agreement with previous studies, we find that the ra- 
dial profile of subhalos within the virial radius in the dis- 
sipationlcss simulations is significantly shallower than the 
dark matter distribution, if the subhalos are selected based 
on the present-day values of the subhalo mass, M°. Com- 
paring two simulations of the same cluster with an order 
of magnitude different resolution, we show that both the 
radial profile of subhalos and DM profiles have converged 
at the scales analyzed here. The difference between matter 
and subhalo distributions is therefore real and is not due to 
numerical overmerging, in agreement with conclusions of 
Dicmand ct al. (2004) . Selection on the present-day value 
of the maximum circular velocity, V^, results in steeper 
radial profile of subhalos. The radial bias with respect to 
the dark matter distribution in this case is weaker than in 
the case of mass-based selection. 

We show that the radial bias of subhalo distribution is 
due to the tidal mass loss experienced by subhalos as they 
orbit in the cluster potential. We find that subhalos in the 
inner regions have on average suffered larger tidal mass loss 
than the halos near the virial radius (see also Gao et al. 
2004b). For instance, the halos within 0.3r v i r on average 
have lost more than 70% of their mass since accretion, 
while subhalos at r > 0.5r v i r on average lose only < 40% 
of their mass. Mass-based selection, therefore, biases the 
subhalo sample to large radii where tidal stripping have 
depleted the number of objects in a given mass range to 
a lesser degree. The mass loss is also accompanied by de- 
crease in the maximum circular velocity, but the decrease 
is slower compared to mass. The slower evolution of the 
maximum circular velocity compared to mass means that 
when it is used to select subhalos the radial bias is smaller. 

One of the obvious omissions in dark matter only simu- 
lations is the lack of dissipation. One can thus ask whether 
the baryon dissipation makes the subhalo more resistant to 
tidal disruption. If so, is the effect strong enough to resolve 
the differences between matter and subhalo distribution? 
Comparing the simulations with and without cooling and 
starformation, we find that the baryon dissipation indeed 
increases the survival of subhalos, especially in the inner 
regions of clusters. However, the effect is relatively small 
because the baryon dissipation makes the core more resis- 
tant to the tidal disruption, but not the entire structure 
of the subhalo. The efficiency of tidal losses for the total 
mass is, therefore, largely unaffected by the baryon dis- 
sipation except in the center of the subhalos. Thus, we 
conclude that this effect alone does not resolve the differ- 
ences between matter and subhalo distributions. 

If the radial bias is due solely to the varying amount 
of mass loss at different radii, the bias can be minimized 
if we use a property of halos unaffected by the evolution. 
In dissipationless simulations such property can be com- 
puted at the actual time of accretion if the evolution of 
a halo is traced in time. For instance, the radial pro- 
files of subhalos selected using mass or maximum circular 
velocity measured at the time of subhalo accretion onto 
cluster are very close to the cluster DM profile. The radial 
distribution of subhalos and inferences from comparisons 



with observed radial distribution of galaxies in clusters 
therefore depend on the selection criteria used to define 
subhalo sample. Consequently, the differences between 
subhalo and galaxy radial distributions do not necessarily 
indicate that dissipationless simulations significantly un- 
derestimate the number of galactic halos in clusters (e.g., 
Diemand ct al. 2004; Gao et al. 2004a). They may simply 
indicate a difference in selection criteria for subhalo and 
galaxy samples. 

Unfortunately, halo properties at the time of accretion 
are not observable in reality. The stellar mass, however, 
is an observable that can behave similarly. Stellar mass 
of a galaxy should change little until it is almost fully de- 
stroyed by tides because stars are located in the centers of 
halos and are tightly bound. The total mass, on the hand, 
can decrease dramatically due to tidal stripping. Using 
the simulations that include cooling and starformation, we 
show that the radial distribution of subhalos selected based 
on the stellar mass is considerably steeper than the that 
of subhalos selected using the total masses in the same 
simulation. 

We show that the profiles of galaxies in the simulations 
are in good agreement with the observed projected distri- 
bution of galaxies for the entire radial range probed by 
the current simulations: 0.1 < R/R 20 o < 2.0. The NFW 
fits to the average profiles of eight simulated clusters at 
z = in the range 0.1 < R/R200 < 1-0 give the best-fit 
concentration parameter of CGai ~ 2 — 3 for galaxies and 
Cdm ~ 10 for the dark matter. In other words, the radial 
profile of galaxies selected using the stellar mass is, on 
average, somewhat more extended than the dark matter 
profile. Note that the values of CQai in our simulations are 
consistent with the observational estimates of CGai ^3 — 4 
(Carlberg et al. 1997; van der Marel et al. 2000; Lin ct al. 
2004). There are indications that in group-size systems 
galaxies have more extended radial profiles than DM (van 
den Bosch et al. 2004; Mathews et al. 2004). The con- 
centration of DM profiles are expected to be in the range 
cdm ~ 5 - 10 (Navarro et al. 1997; Bullock et al. 2001; 
Eke et al. 2001; Tasitsiomi et al. 2004) or larger if the 
profiles are significantly affected by gas cooling (Gncdin 
et al. 2004). We do not find any strong correlation be- 
tween galaxy and DM concentrations in simulated clusters 
(see Figure 10). Thus, our results imply that the projected 
radial profiles of galaxies in clusters may, in general, be a 
poor probe of the underlying dark matter distribution. 

The fact that there is a difference between galaxy and 
DM radial profiles indicates that galaxies are affected by 
tides in the dense environment of the cluster cores. Higher- 
resolution simulations will be needed to test whether the 
amount of the stellar mass loss and efficiency of tidal dis- 
ruption have converged. If the flattening of the inner pro- 
file of galaxies is due to real tidal mass loss or disruption 
experienced by galaxies, we expect to find a large amount 
of stellar debris in the cluster core. Therefore, the detailed 
studies of the properties and origin of the stellar debris and 
intracluster light should provide new insights into the con- 
nection between the evolution of galaxies in clusters and 
the formation of cD galaxies and their stellar envelopes. 
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APPENDIX 
NFW PROFILE FITTING PROCEDURE 

For each halo, we fit the NFW analytic density profiles 
of the form 



p(r) = 
S(r) = 
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where with x = r/r s . The former is the 3D analytic NFW 
profile. The latter is the analytic 2D profile obtained by 
projecting the 3D profile along the linc-of-sight from the 
negative to positive infinity and f(x) is given by (Bartel- 
mann 1996): 
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The NFW profile has two free parameters r s and p s or, 
equivalently, the concentration parameter c v i r = r v ; r /r s 
and the virial radius r v ; r of the halo. These parameters 
are highly degenerate, which makes the fits of the analytic 
profiles to the profiles of the simulated clusters sensitive 
to a number of factors, including the choice of binning, 
the merit function, the range of radii used in the fitting, 
the weights assigned to the data points. As we are only 
interested in the concentration parameters of galaxies and 
dark matter of the simulated clusters, we can improve the 
quality of the fits by removing the degeneracy between 
r v ; r and c. This is possible because we can measure r v i r 
directly in the simulations, leaving concentration as the 
only one free parameter. 

When fitting the profiles of the simulated clusters using 
the one-parameter NFW analytic profile, we first normal- 
ize the radius to the virial radius, y = r/r v i r . For the 



12 




R/R,, 



Fig. A12. — Middle panel: Results of the NFW fits to the pro- 
jected radial distribution E(iJ)/E„; r of dark matter and galaxies 
averaged over the eight simulated clusters at z = and three orthog- 
onal projections. The fits were done using the range [0.1r v ; r , r v ; r ] 
and x 2 merit function. Three lines indicate dark matter (solid) 
and galaxies selected using the thresholds of M* > lO 10 /! -1 *^ 
(dotted) and Aftot > 2 X 10 11 h~ 1 M ( r ) (dashed). Bottom panel: 
fractional deviations of each best fit NFW profile from the actual 
profile. Top panel: local logarithmic slope as a function of ra- 
dius for the three fits. The distribution of both dark matter and 
galaxies is well-described by the NFW profile in the radial range 
0.1 < R/Raoo < 1-0. We find that the radial distribution of galax- 
ies is shallower than that of dark matter, with the concentration 
parameter of 6.3 for dark matter and concentration of 2.0 and 0.87 
for subhalos selected using M* and Mtot, respectively. 



the cluster, beyond which the presence of strong fluctu- 
ations due to cluster substructure significantly biases fits 
of smooth analytic profiles. We checked that the results 
are not sensitive to the variation around the chosen num- 
ber of bins and radial range. For example, decreasing the 
lower/upper bound by a factor of two change the resulting 
best fit concentration of dark matter by only about 10%. 
Note that the insensitivity of the results to the perturba- 
tion around the chosen fit parameters is partly because 
we removed the degeneracy between c and r v i r by measur- 
ing r v i r directly from the simulations. In fact, the same 
changes affect the results of the two-parameters fits much 
more dramatically. Figure A12 shows the best-fit concen- 
tration parameters of the dark matter and galaxy distri- 
bution averaged over eight simulated clusters at z = 0. 
We find that the concentration parameter is 6.3 for dark 
matter, 2.0 for galaxies (selected using M*) and 0.87 for 
subhalos (selected using M to t)- Note that the average con- 
centration parameter of DM obtained here is lower than 
the average value of « 10 quoted in § 7. This is because 
fitting of the projected profiles generally biases the profile 
due to the projection from the large distances outside the 
virial radii, as we did not attempt to subtract the back- 
ground in these two dimensional fits. 



3D profile, we then need to normalize the analytic pro- 
files in units of the mean mass density within the virial ra- 
dius, p vir = 3M v i r /47T7^ ir . Similarly, for the projected (2D) 
profile, we normalize the profile in units of the projected 
surface mass density enclosed within the sphere with the 
virial radius, S v ir = M v i r /irr^ ir . With this normalization, 
the NFW analytic density profiles in Eq. (Al) become 

p{r) _ 1 cZgjc)- 1 ^ 

p vir 3 cy(l + cy) 2 ' 

s vir 2(c V y-i I{cvh 

with only one free parameter, c. Here, g{c) = ln(l + c) — 
c/(l + c). 

For each halo, we obtain the best fit concentration pa- 
rameter c by fitting the analytic profiles given by Eqs. (A3) 
to the 3D and 2D radial profiles of dark matter, galaxies 
and subhalos in the simulations by the standard \ 2 min- 
imization. We use equal-size logarithmic bins and weigh 
each bin by the Poisson noise in the number of dark mat- 
ter particles, galaxies or subhalos (see Tasitsiomi et al. 
2004, for more details). For both dark matter and galax- 
ies/subhalos, we perform fits using the 12 equal-size log- 
arithmic bins between 0.1< r/r V i r <1.0. The fitting re- 
gion is chosen to maximize the statistics of galaxies, since 
the radial profiles of galaxies are much noisier than the 
dark matter profile. The lower radial bound is set by 
the fact that the halo finder fails to find galaxies within 
r/r v i r < 0.1, while the outer bound is the virial radius of 



